function bciCheckFeatures(runs,method)
% bciCheckFeatures(runs,method)
% Plots impact of features 
% INPUT:
%   runs - vector of runs to involve in calculation
%   method - string of comparison method:
%             'tval' - t-values
%             'rsqu' - Rsquare values
%             'svm' - weights of a support vector machine
% OUTPUT: 
%   suggestion of frequencies regarding predefined bands
%
% CR wrote it

global GLOBALbci;
global GLOBALtrainDat;

bci=GLOBALbci;
% here we do not need a featureselection
bci.param.featureSelection='';
[bci,Train]=bciGetTrainDat(GLOBALbci,GLOBALtrainDat{runs});

frequencies = bci.param.freq(bci.init.featOfInterest);
if length(frequencies)~=bci.reshape.reSize(2),
    error('Frequency mapping not unique. bciRefreshParam may help.');
end
if nargin<2,
    method='tval';
end
if nargin<3,
    C=[];
end

if strfind(lower(method),'tval'),
    methodStr = 'Tvalues';
elseif strfind(lower(method),'rsqu'),
    methodStr = 'Rsquare values';
elseif strfind(lower(method),'svm'),
    methodStr = 'SVM weights';
end

xT=1:length(frequencies);
xT=1:ceil(length(xT)/20):length(xT);
xTL={int2str(round(frequencies(xT))')};
cond=bci.eventsToClassify;
nPairs = (length(cond)-1)*length(cond)/2;
meanW = zeros(nPairs,bci.reshape.reSize(2));
maxW = zeros(nPairs,bci.reshape.reSize(2));

W = bciMultiCompareFeatures(bci,Train,method,2);

figure;
% one vs one
k=0;
for k1=1:length(cond)-1,
    for k2=k1+1:length(cond),
        k=k+1;
        meanW(k,:)=mean(abs(W{k}),1);
        maxW(k,:)=max(abs(W{k}),[],1);
        subplot(nPairs,1,k);
        %imagesc(abs(W{k}));
        maxAbsW=max(abs(W{k}(:)));
        imagesc(W{k});
        xlabel('Frequency [Hz]');
        title([methodStr ' ' int2str(cond(k1)) ' vs ' int2str(cond(k2)) ]);
        ylabel('Channel #');
        %colormap copper
        set(gca,'clim',[-maxAbsW maxAbsW]);        
        colorbar
        set(gca,'xtick',xT,'xticklabel',xTL);
    end
end

% determine frequencies in predefined bands
bands=[4 8; 8 10; 10 15; 15 30; 30 40; 40 50; 50 60; 60 80;80 100];
maxIdxBin=zeros(size(maxW));
meanIdxBin=zeros(size(meanW));
meanTthresh=2;
maxTthresh=max([2*ones(size(maxW,1),1),max(maxW,[],2)*0.7],[],2);
for k=1:size(bands,1),
    bandIdx=find(frequencies>bands(k,1)&frequencies<=bands(k,2));
    if ~isempty(bandIdx),
        for c=1:size(maxW,1),
            [m mIdx]=max(maxW(c,bandIdx));
            maxIdxBin(c,bandIdx(mIdx))=maxW(c,bandIdx(mIdx))>maxTthresh(c);
            [m mIdx]=max(meanW(c,bandIdx));
            meanIdxBin(c,bandIdx(mIdx))=meanW(c,bandIdx(mIdx))>meanTthresh;
        end
    end
end
disp('Suggested frequencies:');
informativeFreq=unique(round(frequencies(sum(maxIdxBin,1)+sum(meanIdxBin,1)>0)));
disp(informativeFreq);
